Inherent Rheology of a Granular Fluid in Uniform Shear Flow 
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In contrast to normal fluids, a granular fluid under shear supports a steady state with uniform 
temperature and density since the collisional cooling can compensate locally for viscous heating. It 
is shown that the hydrodynamic description of this steady state is inherently non-Newtonian. As a 
consequence, the Newtonian shear viscosity cannot be determined from experiments or simulation 
of uniform shear flow. For a given degree of inelasticity, the complete nonlinear dependence of 
the shear viscosity on the shear rate requires the analysis of the unsteady hydrodynamic behavior. 
The relationship to the Chapman-Enskog method to derive hydrodynamics is clarified using an 
approximate Grad's solution of the Boltzmann kinetic equation. 
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I. INTRODUCTION 



Consider a fluid between parallel plates in the x, z 
plane separated by a distance L and with relative veloc- 
ity U along the positive x direction. Given appropriate 
boundary conditions, the fluid undergoes simple shear 
with the local velocity fleld given by — ay, Uy = Uz = 
0, where a — U /L is the constant shear rate. The work 
done by the plates on the fluid near the boundaries tends 
to increase the temperature locally due to its viscosity. 
For normal fluids, this is compensated by a heat flux to- 
ward the center with a corresponding temperature gra- 
dient characterizing Couette flow in the steady state. A 
quite different steady state is possible for granular fluids 
where both the temperature and density fields are spa- 
tially uniform, called simple or uniform shear flow (USE) 
P, I2] ■ This is possible because the particle collisions in 
a granular fluid are inelastic and there is a continual 
loss of energy. This collisional cooling can compensate 
locally for the viscous heating so that no heat flux is 
generated. If the temperatures at the walls are not con- 
trolled the fluid autonomously seeks the temperature at 
which this exact balance between collisional cooling and 
viscous heating occurs. Otherwise, for flxed wall tem- 
perature one or the other will dominate leading again to 
Couette flow but with the curvature of the temperature 
field controlled by both mechanisms. This more general 
case has been discussed in detail elsewhere Q including 
USF as a special case. Here, attention will be limited to 
the special features of USF. 
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The steady USF has been extensively studied. Molec- 
ular dynamics and Monte Carlo simulations 0, 0, 0, S 
IE I&3 Jill have been performed to measure the depen- 
dence of the stress tensor on the density and inelastic- 
ity. On the theoretical side, Lun et al. have ob- 
tained the Theological properties of a dense gas for small 
inelasticity, while Jenkins and Richman have used 
a maximum-entropy approximation to solve the Enskog 
equation. This method has been subsequently extended 
to highly inelastic spheres ■ Sela et al. ~\A\ have solved 
the Boltzmann equation to third order in the shear rate, 
finding normal stress differences. Some progress has been 
made by using model kinetic equations for dilute gran- 
ular gases jiy, as well as for dense granular gases [T6| . 
Exact solutions derived in both cases compare quite well 
with Monte Carlo simulations, even for strong dissipa- 
tion. Similar studies for multi-component systems are 
much scarcer, although some recent work has been car- 
ried out [lilillillia. 

Normally, transport properties are defined for general 
states with the temperature, shear rate, and restitution 
coefficient as independent variables. For example, the 
shear viscosity as given by the Chapman-Enskog solu- 
tion to the Boltzmann equation "2(11 has the form for 
shear flow 77 = 77 {T,a,a), and the Newtonian shear vis- 
cosity is obtained from it at zero shear rate, ryo (T, a) = 
77 (T, a = 0, a). However, in the steady state of USE the 
condition for balancing collisional cooling and viscous 
heating implies that the steady state temperature is a 
function of the shear rate a and the coefficient of resti- 
tution a measuring the degree of inelasticity of the par- 
ticles, Ts — Ts{a, a). In this case rj = rj {Ts{a, a), a, a) = 
r]s{a,a) = ri'^{Ts,a). This complication raises two inter- 
esting questions: 

1. For sufficiently small a, is it possible to describe 
USE using Newtonian hydrodynamics? 

2. If ?7s(a, a) is measured in USE, can the Newtonian 
viscosity rjo (T, a) be deduced from it? 
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The objective of this paper is to answer these two ques- 
tions in the general context of hydrodynamics for granu- 
lar fluids. It is shown that the answer is negative in both 
cases, implying that USF is inherently non-Newtonian 
and that the full nonlinear dependence of r]{T,a,a) on 
the shear rate is required, even for small a, at any a ^ 1. 
This result is important because it identifies a fundamen- 
tal and unbridgeable gap between Navier-Stokes hydro- 
dynamics and the hydrodynamics for steady shear flow. 
There is no possibility of using USF to study the Newto- 
nian viscosity (e.g., by molecular dynamics simulation); 
there is no possibility of using the Navier-Stokes equa- 
tions to describe USF or states near it (e.g., the stability 
of USF). This is not an entirely negative result, however, 
because it shows that USF is a rich testing ground for 
the study of rheology. Any measurement will necessarily 
be illustrating some rheological effect. 

These results may seem counter-intuitive at first as 
with many conceptual differences between granular and 
real fluids. However, the steady state considered here 
for the granular fluid does not exist for a real fluid so 
intuition is not reliable in this case. The analysis pre- 
sented here, and illustrated in detail for a dilute gas, is 
important because similar failures of the Navier-Stokes 
level hydrodynamics must be expected for other steady 
states specific to granular fluids. This peculiar limita- 
tion of Newtonian hydrodynamics may apply generically 
when the steady state exists only because of the internal 
cooling mechanism, rather than as a consequence of con- 
trollable external fields or boundary forces. The steady 
state implies the equivalence of spatial heat conduction, 
viscous heating, and coUisional cooling, and a consequent 
relationship between the heat flux, momentum flux, and 
cooling rate (see below). For a given cooling rate (i.e., 
given restitution coefficient) the gradients of the hydro- 
dynamic fields can no longer be controlled independently 
to assure the validity of the Newtonian limit. This phe- 
nomenon is elaborated further in the next section for USF 
based on the macroscopic balance equations for mass, 
energy, and momentum, and idealized boundary condi- 
tions for a general granular fluid. This analysis assumes 
only the existence of a hydrodynamic description, but 
is not limited to any small gradient, small inelasticity, 
or other approximation. A more detailed illustration 
from the Boltzmann kinetic theory is given in Section 
mil An important point of this kinetic theory analysis 
is the determination of the rheological equation of state 
77 (T, a, a) for both steady and unsteady states of USF. In 
the unsteady state the temperature increases (decreases) 
depending on whether the viscous heating of the initial 
state is larger (smaller) than the coUisional cooling. The 
three parameters a, a, T can be controlled independently, 
and there is a small shear rate domain for which the un- 
steady fluid is Newtonian. In the steady state, however, 
these parameters are constrained, T = Ts{a, a), such that 
ri{Ts{a,a),a,a) is never in its Newtonian limit for any 
a < 1. These results and their importance for the study 
of transport in granular fluids are summarized with com- 



ments in the last section. 



II. HYDRODYNAMICS FOR USF 

An idealized granular fluid consists of smooth hard 
spheres {d = 3) or disks {d = 2) of diameter a and 
mass m. The collisions between particles are character- 
ized through a constant coefficient of normal restitution 
a with values < a < 1. The elastic limit corresponds 
to a = 1. The exact macroscopic balance equations for 
mass, energy, and momentum are 

A?^ + 7lV•u = 0, (1) 



(A + C)r+-^(P„V,M, + V-q) =0, (2) 



DtU^ + {mn)-^\/jE,j ^0, (3) 

where n is the density, T is the granular temperature, 
u is the flow velocity, Dt = 9t -|- u • V is the material 
derivative, C is the cooling rate, P is the pressure tensor, 
and q is the heat flux. On physical grounds, C has essen- 
tially the form C « ^{l — a^), where v (x Vt is a mean 
collision frequency and 1 — is the fraction of energy 
lost in each inelastic collision. In order to get a closed 
set of hydrodynamic equations from (Q)^®? Ci ^-nd q 
must be specified as functionals of the fields n, T, and 
u. However, some interesting results can be obtained for 
the purposes here even at this exact level. We empha- 
size that the analysis of this section and its conclusion do 
not depend on any specific form for the rheological equa- 
tion of state, only its scaling with respect to the available 
hydrodynamic fields (dimensional analysis). 

An idealized macroscopic state of USF is characterized 
by the forms 0| 

n{r,t)=n, T{r,t)^T{t), ^ ay, Uy^u^^O. 

(4) 

The presumed geometry is that described above. More 
precisely, this linear velocity profile assumes no boundary 
layer near the walls and is possible for special periodic 
boundary conditions in the local Lagrangian frame '22l|. 
Since the density is a constant, it plays no significant 
role in the following and the dependence of properties on 
it will be suppressed. From symmetry, the cooling rate, 
heat flux, and pressure tensor must have the forms in the 
hydrodynamic state 

C^C{Tit),a,a), q = 0, V./^y - 0. (5) 

The balance equations (QJ and (|3J| are satisfied exactly 
with these choices, while Eq. becomes 

dtT{t) = ~^aP,y {T{t),a, a) - ( {T{t),a, a) T{t). (6) 
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The functional forms for the scalars C (T, a, a) and 
Pxy (T,a,a) must be determmed from a more micro- 
scopic basis. An example from kinetic theory is provided 
in the next section. Then, Eq. JHJl provides a closed hy- 
drodynamic equation to determine the temperature T{t) 
for given initial condition r(0) and specification of the 
constants a and a. The first term on the right side is 
positive and represents viscous heating. To make this 
more explicit it is usual to introduce the nonlinear shear 
viscosity rj (T, a, a) by 



Pxy (T, a, a) = -77 (T, a, a) a. 



(7) 



A Newtonian fluid is that for which the viscosity and 
cooling rate becomes independent of the shear rate, 

T] (T, a, a) ^7] (T, a = 0, a) = ?7o (T, a) , 
C(T,a,a)^C(T,a = 0,a) = Co(T,a). (8) 

The corresponding hydrodynamic equation for a Newto- 
nian fluid becomes 



dtT{t) = -Co (Tit), a) T{t) 
2 



-Ha" 



dn 



7^oiT{t),a)~C2{T{t),a)Tit) 



.(9) 



Here, for consistency, the cooling rate also has been ex- 
panded to second order in the shear rate. 



C {T{t), a, a) -> Co {T{t),a) + a^C2 {T{t),a) . 



(10) 



The Newtonian fluid equation ^ appears quite useful 
as a means to measure its viscosity. First, the cooling 
rate Co {T{t), a) is determined by measuring the time de- 
pendence of T{t) at different values of a at zero shear 
rate. Next, the same measurement is performed as a 
function of the shear rate to determine the combination 
(2/dn) rjo (T(t), a) - (2 a) T{t). In the case (2 = 

(see the next section for an example) this determines the 
Newtonian viscosity. Two questions arise at this point. 
First, can the measurement of the Newtonian viscosity be 
simplified by considering only the steady state? Second, 
under what conditions does the Newtonian limit for the 
steady state apply? To answer these questions we have to 
take into account that the domain of validity of the New- 
tonian description is restricted to small shear rates. To 
be more precise, a dimensionless shear rate must be in- 
troduced. The only other relevant frequency for the fiuid 
is its mean collision frequency i/(T(t)) cx ^/T{t). Con- 
sequently, the relevant dimensionless shear rate, cooling 
rate, and shear viscosity are defined by 

a*{T{t)) ^ C (a ,a)- v : ; 



HT{t)y ' ' ' ' u{T{t)) ' 

r,*(a*,a) = ^^M^(T(t),a,a), (11) 
ni (t) 

where C* (a*, a) and 77* (a*, a) can depend on time and 
temperature only through a*{T{t)), in order to be dimen- 
sionless. Note that the reduced shear rate a* is actually 



a measure of the granular temperature since T (x a* 
In the steady state, Eq. ^ becomes 



*2 



dC* {a*, a) 



(12) 



27/* [a*, a)' 

The solution to this equation gives a* (a), that depends 
only on a and thus cannot be made small by controlling 
the shear rate a. In general, al{a) <C 1 only in the 
quasielastic limit 1 — a ^ 1, 



of (a) 



a-o (a) 



l-a<l, (13) 



where we have taken into account that Coi^) oc ^ ~ 
c? ~ 2(1 — Of) in the quasielastic limit. In this case 
the dimensionless steady state shear viscosity becomes 
77* (a) = 77* (a* (a), a) if {ao(a\ a) . Since, according 
to Eq. ifT^ . Og is small, the steady state shear viscosity 
is expected to differ from the dimensionless Newtonian 
viscosity 770(0) = '7*(a* = 0, a) by a (super-Burnett) 
term of order a?;^: 



(14) 



On the other hand, by expanding the Newtonian shear 
viscosity 77q(q!) around its value in the elastic case, we 
have 

= 1) ~ 1 - (15) 

for small inelasticities. Therefore, Eqs. (|14|l and H15|l 
show that in the quasielastic limit the steady state shear 
viscosity 77* (a) differs from the Newtonian viscosity 77Q (a) 
as much as the latter differs from that of the elastic gas. 
This shows that a measurement of the steady state tem- 
perature to determine a* (a) does not allow determina- 
tion of the Newtonian viscosity, even for 1 — a ^ 1. 
The a dependence of 770(0;) = 77* (a* — 0,q;) cannot be 
isolated from the a dependence of 77* [a\(a),a). Both 
viscosities coincide only in the trivial case a — \. 

In summary, for any chosen a, a, r(0) the steady state 
viscosity is r\* (a* (a), a) and the system samples a non- 
Newtonian value from the general rheological equation of 
state 77* (a* , a) that is always different from the Newto- 
nian value. 



III. DESCRIPTION FROM THE BOLTZMANN 
KINETIC THEORY 

To explore this phenomenon in more detail it is nec- 
essary to calculate the rheological equation of state 
77* (a*, a). This requires a more microscopic analysis such 
as kinetic theory. At low density the granular Boltzmann 
equation provides the appropriate starting point. 



A. Boltzmann equation and Newtonian viscosity 

The Boltzmann kinetic equation determines the prob- 
ability density /(r,v,t) for a particle to have position r 
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and velocity v at time t in a low density gas. It has the 
form 



^+v.V)/(r,v,i) = J[r,v|/(i)]. (16) 



The right side describes the effects of inelastic pair colli- 
sions. The detailed form of the Boltzmann collision oper- 
ator, J, is not required here beyond noting the properties 
necessary for the macroscopic balance equations 



civ I V ^\ J[r,v|/(i)] 

im (v — u) 




(17) 



Here the density n, the temperature T, and the macro- 
scopic flow velocity u are defined in terms of /(r, v, t) by 



Here Q (a) = (^o /^^ is the dimensionless cooling rate in the 
homogeneous cooling state, while /? (a) is a dimensionless 
function of the restitution coefficient given in terms of the 
solution to the linearized Boltzmann equation. Explicit 
results for Q and /3 can be obtained by considering the 
leading terms in a Sonine polynomial expansion. In that 
approximation, 



Co(«) 



d + 2 
M 



(l-a^). 



(3 (a) 



1 + a 



2d 



-(!-") 



(25) 



(26) 



Uniform shear flow 



n(r, t) 
n{r, t)u{r, t) 
ln{r,t)T{r,t) 



dv I V ^1 /(r,v,t). 

im (v — u)^ 



(18) 

The two zeros on the right side of Eq. H17() correspond to 
conservation of mass and momentum, while the last term 
results from non-conservation of energy. The properties 
(|17|l lead to the macroscopic balance equations CO"®; 
with the following microscopic expressions for the pres- 
sure tensor, the heat ffux, and the cooling rate: 



Py(r,i) =m j dwV^Vjf{r,x,t), 
q(r,f) = dvF2v/(r,v,i), 



(19) 



(20) 



It becomes prohibitively difficult to go beyond Navier- 
Stokes order to get the exact full dependence of 77* (a* , a) 
on a* from the Boltzmann equation. On the other hand, 
a good estimate is provided by considering the leading So- 
nine approximation to the distribution function (Grad's 
method) ^2M- 

In the special case of USE the solution to the Boltz- 
mann kinetic equation is spatially uniform when ex- 
pressed in terms of the velocity relative to the local flow, 
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ay, Vy 



Vy, 



mann equation H16|l becomes [21 



Consequently, the Boltz- 



dt - aVy 



d 



/(V,i) = J[V|/(t)]. 



(27) 



Multiplying both sides of Eq. (|27|l by niViVj and inte- 
grating over velocity, we get 



C(r,t) = (l-«2)- 



mn 2 a" 



Adr{^)n{r,t)T{r,t) 
^ j d^r j dvi|v-vi|3/(r,v,i)/(r,vi,i). (21) 

In Eqs. ((T^ and ((JOJ, V = v - u(r,<) is the pecuhar 
velocity. 

It is straightforward to determine the Navier- 
Stokes shear viscosity coefficient ryo (T, a) by using the 
Chapman-Enskog method [IJ. The result is [20ll2ll25t 



nT 
KT) 



where 



{d + 2)T{d/2) \m 
is an effective collision frequency, and 



1/2 



Vo (a) 



/3(«) + ^Co(«) 



(22) 



(23) 



(24) 



dtP.j +a{S,,Py, +Sj,P,y) = m / d^rV,VJJ[Y\J] 



= -A,; 



(28) 



The exact expression of the collision integral is not 
known, even in the elastic case. However, a good estimate 
can be expected by using Grad's approximation 



/(V) ^ /o(V) 
where 



1 



P 

2r \nT 



(29) 



/o(V) = n(m/27rT)''/2 exp(-my V2r) (30) 

is the local equilibrium distribution function. When Eq. 
(|29|l is inserted into the definition of A^ and terms non- 
linear in Pij/nT — 5ij are neglected, one gets |2^ 

A,,=ly[|3iP,J~S,,) + QP^J], (31) 

where Q and /3 are given by Eqs. and respec- 
tively. It is worth noting that Eq. H31|) coincides with the 
one obtained from a simple kinetic model [Tsl |27| . 
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The three relevant independent equations from 128() 
and H31|l are 



2a 

dtp + CqP + -jP^v = 0' 



(32) 



dtPccy + iPi^ + Co) P.y + aPyy = 0, (33) 



As anticipated above, it is independent of the initial tem- 
perature and shear rate a. Appendix fXl shows that the 
steady state solution is indeed a (hnearly) sta- 

ble solution. 

The second equality in H41I) allows one to identify the 

steady state shear viscosity as 



p. 



xy 



^P{a) [P {a) + Q{a)Y 



(43) 



dtPyy + {pV + Co) Pyy - PVP = 0. 



(34) 



where p — nT — Pu/d is the low density pressure. If we 
define the dimensionless quantities 



pm = 



nT{t) 



, a*{t) = 



v{T{t)y 



Tit) = / dt'v{T(t')), 



then Eqs. lI221)-(IS3|l become 



2dr\na* =Q + -;P:ya\ 



(35) 



(36) 



This gives the explicit form for the steady state shear vis- 
cosity, showing it is also independent of the shear rate a 
and the initial temperature. Furthermore, its dependence 
on a is qualitatively different from that of Eq. (|24|l for a 
Newtonian fluid. This is illustrated in Fig. ^ which also 
shows that the analytical results compare favorably well 
with simulation data obtained from the DSMC method 
|28. 29] . Clearly, there is no relationship of the steady 
state shear viscosity to the Newtonian viscosity at any 
value of a. While the Newtonian viscosity is ?7o(Qf) > 1, 
the steady state viscosity is ry* (a) < 1. In the quasielastic 
limit, Eqs. |(23I), ESI, and ||^ yield 



1 



M 



Ad 



-(1-a), 



(44) 



d P* 



^ Pyy Pxy P 



-P* a 



(37) 



af(a)^af(a) = ^(l-a), (45) 



yy 



-P* n" 



(3. 



(38) 



The variable r is the dimensionless time measured as the 
average collision number. The solution to these nonlinear 
equations gives a*, P*y, and P*y as functions of r for a 
given value of a. The rheological equation of state is then 
obtained from 



(39) 



(46) 



Equations (|44|I - H46|I confirm the qualitative arguments 
behind Eqs. ((T^ - ((T^ . Eliminating a between Eqs. (03)) 
and one has 



in the quasielastic limit. 



10 ,2. . 

a, (a) 



d{d + 2) 



(47) 



1. Steady state solution 



2. Unsteady hydrodynamic solution 



Consider first the steady state solution. Equation H36|) 
gives 



p* 



while Eqs. ^ and (jSHJ give 



(40) 



P;y,.=f3if3+Cor\ Ply,s=-<P{P + Qo) (41) 

The value of a* in the steady state is obtained from Eqs. 
and EU, 



«:(a) = W^||^[/3H+Co(«)]- (42) 



The more general solution to Eqs. (|36|I - H38|I corre- 
sponding to hydrodynamics is that for which all time 
dependence occurs through the hydrodynamic fields (the 
Chapman-Enskog "normal" solution). The time depen- 
dence is due only to the temperature for USF which oc- 
curs in the dimensionless forms above through their de- 
pendence on a*{t) = a/v{T{t)), 



dP* 
■' oa* 

Then, using Eq. (jSHl, Eqs. and become 



dP: 



xy 



-2P* 



- P* n* 

i^xy"- 



do 



Co + dPiy°'* 



(48) 



(49) 
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a. 





FIG. 1: Plot of the dimensionless Newtonian shear viscosity 
rio{a) (dashed Hne) and the dimensionless steady state shear 
viscosity ??s(Q!) (solid line) for d = 3, as given by Eqs. (I24|l 
and 1431 1. respectively. Filled circles gg] and open circles |29ll 
represent simulation data obtained from the numerical solu- 
tion of the Boltzmann equation by the DSMC method. 



FIG. 2: Plot of ?7*(a*,a) as a function of a* for d — 3 and 
a = 0.9 (solid line), a = 0.5 (dashed line), and a — (dotted 
line). The thick dashed line is the locus of points {al,ris), 
which are parametrically found from Eqs. 14211 and 14311 . It 
intercepts the curves representing rj* (a* , a) at the steady state 
values (indicated by circles). Note that the curve {a1,,rig) 
ends at the point corresponding to a = (represented by a 
filled circle). 



da* a*{Co+lP*ya*) 

This is a set of two coupled nonlinear differential equa- 
tions that must be solved with the appropriate boundary 
conditions to get the hydrodynamic solution. There is a 
singular point corresponding to the steady state solution 
(|41|l and (|42|l . in which case the numerators and denom- 
inators of Eqs. and (|Sni) vanish. 

The unsteady hydrodynamic solutions ri*{a*,a) ~ 
—P^y/a* are illustrated in Fig.|21for three different val- 
ues of a. Also shown are the special values of the steady 
state shear rate a* (a) and shear viscosity 77* (a) for each 
curve. For each value of a the point (a*,?/*) splits the 
curve 77* (a*) into two physically different branches, one 
for a* < a* and another one for a* > a*. Suppose 
that for a given shear rate a and initial temperature 
T(0) the value of a* is less than that for the steady 
state, a*(0) = a/v{T{0)) < a*. The cooling domi- 
nates viscous heating in this case and the temperature 
decreases, leading to a larger value of a*(T). In this 
way, the system evolves according to the hydrodynamic 
equations along the curve until the steady state value 
a* is attained, i.e., a*(T) a* as r ^ 00. A para- 
metric plot of 77* (r) versus a* (r) gives the branch of the 
curve r]*{a*) corresponding to a* < a*. Analogously, a 
heating process occurs in the opposite case of an initial 
value for a* greater than that for the steady state un- 
til limT-^oo a*(T) = a* again. This provides the branch 
corresponding to a* > a*. As an example. Fig. O shows 
the time evolution of a*(r) oc 1/-\/T(t), —P*y{T), and 
P* (r) for a = 0.5 and two different initial conditions: 
(i) a*(0) = 0.1, P4(0) = -0.1, P;^(0) = 1, and (ii) 
a*(0) = 2, P* (0) = -0.5, P* (0) = 0.4. In both 



cases, after about 20-30 collisions per particle, the sys- 
tem reaches a common steady state with a* = 0.812, 
-P* = 0.577, and P* , = 0.667. We have checked 
that the same asymptotic state is achieved when starting 
from different initial conditions. 

In order to get the branch of the hydrodynamic solu- 
tion corresponding to a* < a* one must apply the physi- 
cal boundary condition 

lim 7]* {a* , a) = (a) , lim P* (a* , a) = 1 . (51) 

a*— »0 a*— >0 

In practice, the branch a* < a* for each curve of Fig. [21 
has been obtained by solving numerically the coupled set 
(|49|l and (|50() with the initial condition 

P4 = -r;oa*(0), P,; = l (52) 

with a*(0) = 10"^ 

As said before, the branch corresponding to a* > a* 
is inaccessible starting from an initial value a*(0) < a*. 
Therefore, in order to obtain 77* (a*) for a* > a*, one 
must take the boundary condition at the point at infinity. 
For that case, simple physical intuition is not enough to 
determine the appropriate boundary conditions, so one 
must resort to a detailed asymptotic analysis of Eqs. I|49|) 
and H50(l . Such an analysis is carried out in Appendix IbI 
where we find that that -P*y - a*^^/^ P*y - a*"^/^ 
for asymptotically large a* . More specifically, if the shear 
rate a* is much larger than 1, then 
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FIG. 3: Time evolution of the reduced shear rate a*(T), the 
reduced shear stress —P*y(T), and the reduced normal stress 
Pyy(r) for d = 3, a = 0.5, and two different initial conditions: 

(i) a*(0) = 0.1, P^y{Q) = -0.1, Pyy{Q) = 1 (solid lines) , and 

(ii) a*(0) = 2, P*y{Q) = -0.5, Pyy(Q) = 0.4 (dashed lines). 
The horizontal dotted lines represent the steady state values 
a* = 0.812, 



^xy,s — 0.577, 



and Pyy^, = 0.667. 



The branch a* > a* in Fig. |21 has been obtained from 
the numerical solution of the coupled set H49|l and JsT 
with the initial conditions H53|) with a*(0) = 10. 



IV. DISCUSSION 

The primary observation here has been that granular 
fluids admit hydrodynamic steady states that are inher- 
ently beyond the scope of the Navier-Stokes or Newto- 
nian hydrodynamic equations. The reason for this is the 



existence of an internal mechanism, collisional cooling, 
that sets the scale of the spatial gradients in the steady 
state. For normal fluids, this scale is set by external 
sources (boundary conditions, driving forces) that can be 
controlled to admit the conditions required for Navier- 
Stokes hydrodynamics. In contrast, collisional cooling is 
flxed by the mechanical properties of the particles mak- 
ing up the fluid. An example is a sheared fluid where 
the work done at the boundaries is balanced by a com- 
bination of collisional cooling and internal heat flux. To 
illustrate and emphasize these effects in more detail and 
quantitatively, we have considered the idealized state of 
uniform shear flow (USF), for which no heat flux occurs. 
The Lees-Edwards boundary conditions generating USF 
cannot be applied in experimental conditions, although 
they can be easily implemented in computer simulations. 
However, even for the laboratory conditions of Couette 
flow the cooling rate will still be a function of the hy- 
drodynamic gradients leading to similar non-Newtonian 
behavior. 

In spite of the extensive prior work on USF for granular 
fluids, the observation about its inherent non-Newtonian 
character is apparently new although implicit in results 
obtained from various models (e.g., the observation that 
a*^ ~ 1 — a). This is significant because molecular dy- 
namics simulation of steady USF has been used for real 
fluids to measure the Newtonian shear viscosity. The 
analysis here shows that this is not possible for granular 
fluids. More generally, it provides a caution regarding 
the simulation of other steady states to study Navier- 
Stokes hydrodynamics when the gradients are strongly 
correlated to the collisional cooling. On a more posi- 
tive note, these results show that USF is an ideal testing 
ground for the study of rheology since any choice of the 
shear rate and a will provide non-Newtonian effects. It 
is one of the fascinating features of granular fluids that 
phenomena associated with complex fluids are more eas- 
ily accessible than for simple atomic fluids [sOJ. 

The macroscopic state of steady uniform shear flow is 
simple in the sense that it is spatially uniform in the local 
Lagrangian frame, and the only hydrodynamic gradient 
is that of the velocity field characterized by the scalar 
shear rate a. Non-Newtonian effects occur through the 
dependence of the shear viscosity on the shear rate. The 
condition for a steady state is an exact balance between 
viscous heating and collisional cooling. It is shown in 
Section ^] that this implies a value of the appropriate 
dimensionless shear rate that is fixed by the restitution 
coefficient. Hence, it is not possible to make this hydro- 
dynamic gradient small, as required for the Newtonian 
limit, by any initial control of a or the initial tempera- 
ture. Further quantitative illustration of this is given in 
Section ITTTl where the shear rate dependence of the shear 
viscosity is determined from an approximate Grad's solu- 
tion of the Boltzmann equation. The temporal approach 
to the steady state via cooling or heating proceeds along 
a rheological equation of state shown in Fig. [3 but the 
final steady state is independent of the initial conditions 
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and lies in the non-Newtonian domain for any value of 
the restitution coefficient a ^ 1. One consequence is that 
experimental or simulation measurements in steady USF 
provide no information about Navier-Stokes transport. 

This phenomenon of peculiar steady states for granular 
fluids extends to other states as well. The steady state 
macroscopic balance equations are 

V • (nu) = 0, (54) 



Equations lf5S|l and clearly show that A* (a) ^ Aq (a) . 
In the three dimensional case, for instance, the coefficient 
of 1 — a in A* (a) is about twenty times larger than that 
of AQ(a). Therefore, the Navier-Stokes transport coeffi- 
cients do not describe correctly the heat flux in the steady 
state, except again in the trivial case of elastic collisions. 

Acknowledgments 



-u- VlnT 



dnT 



{P^JVJU^ + V • q) = C, 



u 



Vui + (mn) 



(55) 



(56) 



In the absence of collisional cooling the spatial gradients 
can be sustained only by boundary conditions or external 
forces. However, for finite C spatial gradients can exist 
that are not controlled by such external sources and their 
size may preclude the validity of simple Navier-Stokes 
hydrodynamics. Each physical state should be checked 
for this possibility. 

As another simple example, consider the steady state 
of a granular gas enclosed between two parallel plates at 
rest and maintained at the same temperature. For elastic 
particles such a steady state is trivially that of equilib- 
rium at the wall temperature. Nevertheless, for a gran- 
ular gas the collisional dissipation induces a heat flux, 
cf. Eq. (|55|l . and hence a thermal gradient, along with a 
density gradient, exists [^[s^illj- In the quasielastic 
limit, the kinetic model of the Boltzmann equation pro- 
posed in Ref. 527| admits a solution characterized by a 
constant pressure p = nT and a heat flux given by 



d + 2 p 
2 miy 



A:(a)VT, 



(57) 



where 



N . 29d2 + 116d-28,, 

(58) 

On the other hand, the expression for the heat flux in 
the Navier-Stokes order given by the Chapman-Enskog 
expansion is 20, 24, 25] 



q = -Ao(a)VT-^o(a)Vn 

d + P T^, 
= 7, Ao(a)VT, 



(59) 



where Aq is the thermal conductivity, hq is a transport 
coefficient with no analog in the elastic case, and in the 
last step we have taken into account that Vp = 0, so that 
Vn = —{n/T)VT. The expression for the dimensionless 
effective thermal conductivity AQ(a) in the quasielastic 
limit is 



AS (a) 



1-1- 



3d 



-{1-a). 



(60) 
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APPENDIX A: LINEAR STABILITY ANALYSIS 
OF THE STEADY STATE SOLUTION 

In this Appendix we carry out a linear stability anal- 
ysis of the steady state solution H41() - (|42|l of the set of 
evolution equations H36|l - H38() . First, we write 

a*{r) = a*+(5a*(T), P*yiT) = P*y^, + SP*y{T), 
Pyyi^) = PLs + SP;yir). (Al) 

Substituting IjAip into Eqs. H36|) - (|38|l and neglecting non- 
linear terms, one gets 



da* 
SP*y 
<^Pvv 



5a* 
SP* 
a*JP* 



(A2) 



-Co(/3 + Co)V2/3 
P + 2Co 1 



where L is the matrix 
Co/2 

PQ/iP + Q) 

(A3) 

In Eq. (EU use has been made of the explicit expres- 
sions (I41ll - H42|l for the steady state solution. The time 
evolution of the deviations from the steady solution is 
governed by the eigenvalues of L. If the real parts of 
those eigenvalues are positive the steady solution is lin- 
early stable, while it is unstable otherwise. The solution 
of the characteristic equation det(Ly — iSij) = yields a 
real eigenvalue ii and a pair of complex conjugate eigen- 
values £2, £3- We have verified that ii and the real parts 
of (2.3 are positive definite for all values of the coeffi- 
cient of restitution a < 1, £1 being smaller than Re £2,3- 
Consequently, the steady state solution is stable, and the 
characteristic relaxation time (measured by the number 
of collisions) is As an illustration. Fig. ^ shows the 
a-dependence of £1 and Re £2,3 for d — 3. It is inter- 
esting to note that £1 ^ in the elastic limit a ^ 1. 
This is a consequence of the fact that there is no steady 
solution at a = 1, namely a*(T) decays algebraically as 
a*(r) w (2r/d)-i/2. 
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(2/? -62^)^20*""^ =2/?. (B5) 
Equation (|B4p implies that 62 = 6i + 1 > 2, while Eq. 
(|b3|) yields 62 = and A2 = 1. Therefore, Eqs. (|B4|| and 
(|B5(I are mutually inconsistent. The next possibility is 
61 = 1. This gives 



2/3 -Co +1^1 ]Aia*-=2A2a' 



■.-b2 



(B6) 



2p - 02C0 ^ ^ — ^1 



Aaa*"*"' = 2/3. (B7) 



FIG. 4: Plot of the eigenvalue £1 and of the real part of 
£2,3 as a function of the coefficient of restitution in the three- 
dimensional case. 



According to Eqs. (|B6|l and l|B7|) . &2 = 2 and 62 = 0, 
respectively, what is again inconsistent. 

Finally, we consider the case &i < 1, so that Eqs. I)B2|I 
and IIB3I) reduce to 



The above stability analysis is restricted to small de- 
viations from the steady state solution. The proof that 
the steady state solution is an attractor as t — > cx) for 
any uniform (in the Lagrangian frame) initial condition 
seems to be quite difficult, so that one has to resort to 
numerical evidence. As was illustrated in Fig. |31 the nu- 
merical results confirm that the steady state is indeed an 
attractor. 



.2 *-2bi 



A2a , 



These two equations are consistent provided that 



(B8) 



(B9) 



APPENDIX B: HYDRODYNAMIC SOLUTION 
FOR LARGE SHEAR RATES 



h ^ h ^ 



(BIO) 



In this Appendix we will get the asymptotic behav- 
ior of the solution to Eqs. (|49|l and (|50|l for large values 
of the reduced shear rate a*. On physical grounds, the 
asymptotic solution must have the general form 



xy 



~Aia* 



yy 



A.a* 



(Bl) 



where Ai > 0, A2 > 0, and &2 > 0. Inserting (|B1|) into 
Eqs. (03 and ifSUI) . one has 



Ai 



J A2 



V 56 y V 64 

Therefore, for large reduced shear rates, 



p* 



(Bll) 



(B12) 



2(61 



-Aia 



■.l-bi 



Aia*'^'^' = 2A2a*'''' 
(B2) 



yy y 64 



(B13) 



A-,a 



, —bo 



2(3. (B3) 



In order to carry out a balance analysis of the leading 
terms in (|B2|I and (|B3|) . it is convenient to consider sep- 
arately the cases 5i > 1, 61 = 1, and 61 < 1. If 61 > 1, 
then a*^'''- 



(B4) 



and Eqs. l|B2|) and l|B3|) become 



(2/3-6iCo*)Aia 



-1-61 



2yl2a^ 



■-bo 



It is interesting to remark that \P*y\ both for 
small and large shear rates. Thus, \P^y\ exhibits a non- 
monotonic dependence on a* characterized by a maxi- 



duced shear rate. For instance, in the three-dimensional 
case we have found (Omaxj l^^ylmax) — (0.64,0.634), 
(1.22,0.595), and (1.81,0.579) for a = 0, 0.5, and 0.9, 
respectively. A similar non-monotonic behavior of \P*y\ 
occurs in the elastic case as well |21j|. 
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